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ABSTRACT 


This thesis presents the most complete climatology of radio wave (VHF and 
higher frequencies) tropospheric ducting within Turkey and the surrounding areas 
ever published. This ducting can result in greatly extended propagation ranges. 
Atmospheric data from radiosonde upper air observations were collected for a 
five-year period at stations in Adana, Ankara, Diyarbakir, Isparta, Istanbul, Izmir, 
Athens, Bucharest, Crete Island, and Tuapse. The upper level and surface 
ducting yearly average frequencies for these stations were 59.46%, 34.22%, 
37.53%, 38.15%, 50.67%, 45.69%, 64.69%, 30.24%, 65.59%, and 21.05%, 
respectively. Ducts occurred mostly in summer. These frequencies were much 
higher than previously published climatologies. However, many of these ducts 
were too shallow to be operationally significant. Atmospheric pressure was 
correlated with duct heights for some stations, but sometimes the correlation was 
positive while other times it was negative. In most cases, no correlation was 
found between duct height and wind speed or direction. Ducts tended to be 
more common in coastal regions. Regions away from the coast usually have 
topography limits the effects of the ducts away from the coast, and these inland 
regions tended to have less ducting occurrence Therefore, these results are most 
applicable to the coastal regions. 
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I. INTRODUCTION 


The purpose of this study is to find ducting conditions and to quantify 
patterns and trends of them within Turkey and surrounding areas, and to 
understand how ducting can be effective on the performance of communication 
systems. 

Electronic Warfare (EW) plays an important role in the overall 
scheme of military strategy. Information superiority plays a vital role 
in the operation area. A key role of EW is to be able to exploit 
electromagnetic (EM) missions from a hostile source, turn that 
information into intelligence on an adversary’s order of battle, 
intentions, and capabilities, and to render communications and 
weapons systems useless, all while protecting their own system. 
(Cairns-McFeeters, 1992) 

The key matter is how the systems will be used effectively to send and 
receive signals. Many things can affect the electromagnetic radiation in the 
earth’s atmosphere. The ducting is one of the reasons for the causes of 
anomalous propagation. Through ducting, radio frequency detection ranges can 
be extended. Understanding ducting enables one to take advantage of the 
environment in EW. 

Anomalous propagation occurs frequently and therefore must be 
factored into any warfare commander’s decision-making process. 
Weather forecasters supporting military operations should 
understand the effects of atmospheric variables on radio and radar 
performance. They should collect, analyze, and evaluate available 
data to describe the existing atmospheric propagation condition and 
to understand how it will change over time. They should provide 
radar operators and operational commanders with this information 
on a timely basis. (Ford, 2005) 

Determination of refractive conditions can help to predict the following: 

• Flow far away enemy sensors will be able to detect friendly assets 
(missiles, aircraft, ships, boats). 
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• How far away friendly sensors will be able to detect enemy assets, 
and where reconnaissance platforms need to be placed to exploit 
the atmosphere and extend friendly detection capability, 

• What altitudes provides aircraft with the greatest weapons stand-off 
range, and 

• From what direction friendly assets should approach the enemy to 
give the greatest detection advantage (Ford, 2005). 

This thesis includes an overview of EM propagation through the 
atmosphere, discussion of environmental factors that affect EM propagation, and 
finds ducting conditions for ten different locations. It concerns radio waves at 
VHF and higher frequencies and the space wave mode of transmission. 

Atmospheric data was collected from ten different radiosonde stations. 
The stations examined are located in Adana, Ankara, Diyarbakir, Isparta, 
Istanbul, Izmir, Athens, Bucharest, Crete Island, and Tuapse. To find refractivity 
profiles for each radiosonde, a MATLAB program was written. This program 
reads radiosondes, calculates modified refractivity, finds the first three ducts in 
the profile if they exist, and finally does some statistics to evaluate the ducting 
effect on EM propagation. 

This thesis adds power to operators and commanders’ prediction ability. 
Knowing how to use the ducting climatology can extend commanders and 
operators’ line of sight. Communication is like breathe: you have never known its 
importance until it is gone. 
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II. ELECTROMAGNETIC WAVES AND ATMOSPHERIC 

PROPAGATION 


A. ELECTROMAGNETIC WAVE 

Electrons have negative charges and protons have positive charges. If a 
particle is charged, it produces an electric field in all directions. While particles of 
the same charge repel, oppositely charged particles attract each other. If the 
particle moves, the electric field changes. The electrical potential difference 
between two locations is called voltage. 

If a charged particles moves, it creates magnetic field. If an electron 
moves through in a wire, the magnetic field encircles the wire. Its direction 
depends on the direction of electron movement. If a charged particle moves, it 
produces both an electric field (because the particle is charged) and a magnetic 
field (because the particle is moving). An accelerating charged particle creates 
an Electromagnetic (EM) wave. 

Figure 2 shows the creation of EM waves. The magnetic and electric 
fields of an electromagnetic wave are perpendicular to each other and to the 
direction of the wave. The motion of the particle is changing, the electric field is 
changing and the magnetic field is changing. The changing electric field creates 
a new magnetic field, and the changing magnetic field produces a new electric 
field. The collapsing and regeneration of the electric and magnetic fields is what 
allows EM radiation to propagate (Guest, 2005). 


z 



Figure 1. Electromagnetic (EM) Wave Generated by an Accelerating Electron 

[From (Guest, 2005)] 
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B. ATMOSPHERIC PROPAGATION 


Electromagnetic wave propagation is affected by the properties of the 
medium located between the transmitting antenna and the receiving antenna. In 
most situations radio waves in free space travel in straight lines (an exception 
being gravitational bending as predicted by General Relativity theory), while radio 
waves traveling within the earth's atmosphere are affected by varying conditions. 
The influence exerted on radio waves by the earth's atmosphere adds many new 
factors to complicate what at first seems to be a relatively simple problem. 
These complications are due to a lack of uniformity within the earth's 
atmosphere. Atmospheric conditions vary with changes in height, geographical 
location, and even with changes in time (day, night, season, and year). 
Knowledge of the composition of the earth's atmosphere is extremely important 
for understanding wave propagation (Integrated Publishing, 2007). 

Electromagnetic radiation (including radio waves) can be reflected, 
refracted, and diffracted within the atmosphere. To understand radio 
propagation, it is useful to consider different regions of the atmosphere 
surrounding the earth as shown in the Figure 2. 



Ionos])hcre - 


Figure 2. Regions of the atmosphere 


Just above the earth's surface is the troposphere, followed by the middle 
atmosphere (or stratosphere), and finally the ionosphere. The troposphere is 
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characterized by relatively high moisture content and is affected by surface 
topography. Virtually all weather phenomena take place in the troposphere. The 
temperature in this region generally decreases with altitude, clouds form, and 
there may be turbulence due to wind shear or surface heating. Vertical variations 
in temperature, humidity, density, and pressure in the troposphere can have a 
substantial effect on the propagation of radio waves. The stratosphere has much 
more constant conditions and there is little water vapor present. Therefore, the 
stratosphere has relatively little effect on radio propagation. The ionosphere has 
charged particles that cause upward propagating radio waves in the HF and 
lower bands (see Table 1 below) to be refracted back toward the surface. This is 
the most important region of the atmosphere for long distance point-to-point 
communications in the HF range (Davidson, 2003). However, this thesis is 
concerned with tropospheric propagation. 


Table 1. Frequency Ranges and Propagation Characteristics for 

Different Bands 


Bands 

Frequency Ranges 

Propagation Characteristic 

ELF 

0.3 - 3 kHz 

Earth/ionosphere waveguide modes 

VLF 

3 - 30 kHz 

Earth/ionosphere waveguide modes 

LF 

30 - 300 kHz 

Ground waves (Guided by ground) 

MF 

0.3 - 3 MHz 

Ground waves 

HF 

3 - 30 MHz 

Sky waves (Ionospheric “Reflection”) 

VHF 

30 - 300 MHz 

Space Waves 

UHF 

0.3-3 GHz 

Space Waves 

SHF 

3 - 30 GHz 

Space Waves 

EHF 

30 - 300 GHz 

Space Waves 
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How an EM wave behaves in the atmosphere is usually related to its 
frequency. Different frequencies behave differently under the same conditions. 
Table 1 shows frequency ranges and characteristics for different bands. Space 
waves are usually considered to have “line-of-sight” propagation, but refractive 
effects can be very important and increase or decrease the line-of-sight ranges. 

1. Ground Waves 

EM radiation can propagate in various modes of transmission. Ground 
waves are radio waves that travel near the surface of the earth. The ground 
wave is actually composed of two separate component waves: surface waves 
and space waves. A surface wave travels along the surface of the earth. A 
space wave travels over the surface. 

2. Sky Waves 

Sky waves are radio waves that are reflected (actually refracted) back to 
earth from the ionosphere. 



Figure 3. Diagram Showing Ground Waves (Surface and Space) and Sky 

Waves [From (Davidson, 2003)] 
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III. ATMOSPHERIC EFFECTS ON PROPAGATION 


A. ATMOSPHERIC EFFECTS ON RADAR RANGE 


In the atmosphere, electromagnetic waves do not propagate in straight 
paths as they do in free space. They behave differently depending on the 
frequency and the atmospheric effects. This thesis concerns radio waves at VHF 
and higher frequencies and the space wave mode of transmission. In this range, 
the most important atmospheric effects happen in the troposphere, usually in the 
lower part near the earth’s surface. Space waves in the troposphere undergo 
attenuation, scattering, and refraction. 

For radar transmissions, EM energy must propagate toward a target and 
then be reflected back toward the antenna. In free space, radar energy 
propagates according to the simple form of the radar equation (Skolnik, 2001): 


R 


max 




(3.1) 


where: 

Pf = Transmitted power, W 
G = Antenna gain 
Ae = Antenna effective aperture, m^ 
a = Radar cross section of the target, m^ 

Smin = Minimum detectable signal, W. 

In the troposphere, this simple form of the radar equation does not 
adequately predict the range performance of the actual radars. There are 
several reasons for failure of the simple form of the radar equation. This thesis 
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will focus on those situations where and when the simple radar equation and line- 
of-sight propagation do not occur due to effects caused by the earth’s surface 
and atmosphere. 

Propagation effects might extend the radar range or reduce it. It is 
important to consider the earth’s environment when predicting radar 
performance. Scattering (reflection off particles and inhomogenieties), refraction 
(bending), ducting, diffraction, attenuation, external noise and backscatter from 
land, sea and weather clutter, and attenuation in rain and other hydrometeors 
modify free space radar performance and affect the EM propagation (Skolnik, 
2001 ). 

B. REFRACTION 

Unlike free space, the atmosphere is not uniform. Variations in pressure, 
temperature, and most importantly, humidity, cause electromagnetic waves to be 
bent (refraction). A bending toward the surface can extend the radar horizon and 
increases the coverage of the radar, while bending away from the earth will result 
in decreased ranges. Refraction can also introduce an error in the measurement 
of the target elevation angle; the apparent position of a detected target is 
different from its actual position. 
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Figure 4. Tracking Errors for Standard Atmosphere [From (Davidson, 2003)] 

Refraction of the radar waves in the atmosphere is due to the variation of 
the velocity of propagation with altitude. The index of refraction (Equation 3.2) is 
a measure of velocity of propagation in a medium compared to the velocity in 
free space. EM waves travel faster in mediums with relatively lower indices of 
refraction. The index of refraction, n, for a particular medium is determined by 
the magnitude of the phase velocity in space, c, (i.e., the speed of light) divided 
by the phase speed of light in the medium, v. 

n = clv (3-2) 

At the earth’s surface, the value of the index of refraction is around 
1.000315. In order to simplify the use of refraction characteristics it is useful to 
define a quantity called the refractivity, N, which is the index of refraction minus 
one times 10®. This results in easy to use numbers (such as 315 instead of 
1.000315). The refractivity, N, for radio waves depends on atmospheric 
pressure, P (hPa), Temperature, T, (K) and vapor pressure, e, (hPa) according 
to the following formula. 

A^ = (n-l);cl0'' =77.6^-5.6^ + 3.75.xl0'^ (3.3) 
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Note that this formula is valid for all radio frequencies in the VHF or higher 
frequency ranges. The index of refractivity corresponds to N=315 at sea level. 

The vertical gradient of N determines how EM energy will propagate 
through the atmosphere in relation to the horizontal. If dN/dz (N gradient) is 
positive, EM energy is bent upward (subrefraction), and if it is negative, EM 
energy is bent downward. 

To identify trapping conditions and determine the occurrence of ducts, it is 
easy to use a modified index of refraction, M, where z is the height above the 
surface in meters (Equation 3.4): 

M(z) = A^(z) + 0.157Z. (3.4) 

Subicrr.n.'lii>n 



Figure 5. dN/dz and Refraction Categories [From (Davidson, 2003)] 

When dM/dz=0, the ray curvature equals the earth's curvature. This is 
another way of saying that when the N-gradient is -157 units per kilometer, the 
ray has the same curvature as the earth. M will increase with height in the 
standard atmosphere. The importance of M is that for a duct to exist, dM/dZ 
must be negative somewhere in the profile (Davidson, 2003). Figure 5 shows the 
refraction categories according to the N gradient. 

M and N profiles can be compared numerically as shown in Table 2. 
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Table 2. Comparison of N and M Gradients for Refractive Conditions 


Refractive 

condition 

dN/dz 

(N units/km) 

dM/dz 

(M units/km) 

Distance to 

Surface Horizon 

Subrefraction 

0 < N 

157 < M 

Reduced 

Normal 

-79 < N < 0 

78 < M <157 

Standard 

Superrefraction 

-157 < N < -79 

0 < M < 78 

Increased 

Trapping 

N <-157 

M < 0 

Greatly Increased 


A trapping layer is the region where dN/dz < -157 units/km or dM/dz 
< 0. The ray curvature bends downward relative to the earth. A 
trapping layer that does not extend down to the surface is 
sometimes referred to as an elevated layer. A duct is the region 
associated with the trapping layer; duct EM energy is confined and 
channeled between the top and bottom of the duct. The top of the 
duct will always be the top of the trapping layer. However, the 
bottom of the duct can, and often does, extend below the bottom of 
the trapping layer. A duct can have its lower boundary at the 
surface (a "surface based duct") or above the surface (an "elevated 
duct"). (Davidson, 2003) 


Refraction is the bending of radio waves. Unlike reflection, there is no 
sudden direction change in refraction. The EM wave bends as it passes through 
regions with different indices of refraction. The index of refraction of a medium 
changes with the wave’s phase velocity. 


V = 



c 

n 


(3.5) 


As seen from Equation 3.5, the velocity of propagation within any medium 
is a function of the permittivity (dielectric constant) of the medium. The index of 
refraction, n, is equal to^/e^, where is the relative dielectric constant. 
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Figure 6. Snell’s Law 


According to the Snell's Law, when a propagating wave in medium 1 with 
refractive index ni impacts medium 2 with index r \2 at an incidence angle 0i, the 
refraction angle 02 is given by Equation 3.6 (Davidson, 2003): 

n, sin 6*j = sin 0 ^ ■ (3-6) 

Although in the atmosphere n does not change abruptly, Snell’s Law can 
still be used to predict propagation angles resulting from changes in n. From 
Snell’s Law, it is clear that if the 0i is large enough, the critical angle can be 
exceeded and propagation ray paths are almost parallel to the trapping layer. 
The largest possible angle of incidence that still results in a refracted ray is called 
the critical angle; in this case, the refracted ray travels along the boundary 
between the two media. In addition, the vertical position of a target can be wrong 
if there is no correction for atmospheric refraction. 

C. DUCTING 

Normally, in the troposphere, the air is colder at higher elevations. 
However, sometimes the temperature increases with elevation. When this 
happens, it is known as a temperature inversion. Temperature inversions are 
dynamically stable and often prevent the mixing of air by turbulence. This can 

trap moisture at the lower levels and cause a strong decrease in vapor pressure, 
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e, across the inversion. The increase in temperature and more importantly (for 
most situations) the associated decrease in vapor pressure (humidity) can cause 
radio waves to bend toward the earth more sharply than the earth’s curvature, 
and hence the radiation is “trapped” near the surface. This situation associated 
with inversion creates what are known as channels or ducts, where extended 
“over-the-horizon” propagation ranges can occur. If any EM wave reaches and 
enters a duct, this duct acts as a waveguide to send signals long distances. 
Knowing the ducting climatology of a region is important for any Electronic 
Warfare System that transmits or receives radio signals (Davidson, 2003). 
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Figure 7. Temperature Inversion [From 
http://apollo.lsc.vsc.edu/classes/met130] 

The gradient of the modified refractivity with height is the best method to 
see the ducts. Equation 3.3 was about refractivity, N. To identify the trapping 
conditions modified refractivity is used: 

M(z) = /V(z) + 0.157Z = 77.6^^-5.6-^ + 3.75.xl0'-^^ (3.7) 

T(z) T{z) T\z) 
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P = Atmospheric Pressure in milibar 
T = Temperature in Kelvin 
e = Vapor Pressure in milibar 
Z = Height. 



Figure 8. M Relationship for Elevated Trapping Layer and Duct [From 

(Davidson, 2003)] 

Figure 8 shows different types of ducts. The top of the trapping layer is 
the top of the duct. The bottom of the trapping layer can be determined by 
drawing a vertical line from the top of the duct down until it crosses the M profile. 
If it is a surface-based duct, the bottom of the trapping layer is the ground. 

1. Evaporation Ducts 

Evaporation ducts most commonly occur over the ocean and other large 
bodies of water. They are found very close to the surface, in the lowest two 
meters or so. The air right at the surface of a body of water is usually nearly 
saturated with water vapor, which means that the relative humidity, RH, is 100 

14 









percent. Over oceans, salt causes the surface RH to be 98%. Although the 
relative humidity is high over the sea where the air touches, it usually decreases 
rapidly in the several meters just above the sea if the upper air is not saturated 
(i.e., foggy). Because of this rapid decrease in humidity, refractivity also 
decreases rapidly, enough to refract radio waves back toward the surface. Thus, 
a low-lying duct traps the EM waves. Signals can propagate along the surface of 
the water in an evaporative duct. Evaporation ducts can cause the radar ranges 
for targets at or near the sea surface to be considerably greater than the free 
space range; although this effect may be counteracted by the increase in surface 
clutter associated with stronger evaporations ducts (Skolnik, 2001). 



Figure 9. Mirage Over the Water Associated with an Evparation Duct [From 

(Martin, 2007)] 

Figure 9 shows an “inferior” mirage caused by an ocean surface that is 
much warmer than the air. Because warmer water has a relatively higher vapor 
pressure, this situation is also associated with strong evaporations ducts. These 
mirages associated with evaporation ducts are best seen during the daytime, 
though they also occur at night. Evaporative ducts can act as an RF mirror (and 
sometimes an optical mirror) and reflect signals from the top of the duct (Martin, 
2007). 
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2 . 


Surface-based Ducts 


The base of a surface-based duct is at earth’s surface. There are three 
types of surface-based ducts. One of them is the evaporation duct discussed 
previously. The second one is created from an elevated trapping layer, and the 
third one is created from a surface-based trapping layer. 
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Figure 10. Surface Duct for (a) Surface-based Trapping Layer and (b) 
Elevated Trapping Layer [From (Davidson, 2003)] 


As seen in Figure 10, the value of M at the duct top is less than the value 
of M at the surface. 



Figure 11. Trapped Moisture Over Land Due to Radiational Cooling and 
Inversion Formation Near the Surface can Create Surface-based Ducts 
Over Land [From (Martin, 2007)] 
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3. 


Elevated Ducts 


The base of an elevated duct lies above the surface of the earth. This 
type of duct is similar to the surface-based duct. The only difference is that the 
bottom of duct is higher than ground level. This situation is generally caused by 
temperature inversions or decreases in humidity (vapor pressure) with elevation. 



Figure 12. Elevated Inversion with a Base of 620 m is Formed by Subsidence 
in a Fligh-pressure Region and is Often Associated with Elevated Ducts 

[From (Martin, 2007)] 
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IV. SOUNDING DATA SETS AND MATLAB CODES 


A. SOUNDING DATA SETS 

The goal of this thesis is to quantify ducting climatology within Turkey and 
its surrounding areas. All the atmospheric upper air data sets were acquired 
from the Department of Atmospheric Science at the University of Wyoming. Five 
years worth of data were collected to calculate refractivity for ten different 
stations. Table 3 shows these stations and their locations. 


Table 3. Stations and their Locations 


Station Code 

Station Identifier 

Location 

17351 

— 

Adana/Turkey 

17130 

— 

Ankara/Turkey 

17281 

— 

Diyarbakir/Turkey 

17240 

— 

Isparta/Turkey 

17062 

— 

Istanbul/Turkey 

17220 

— 

Izmir/Turkey 

16716 

LGAT 

Athens/Greece 

15420 

LRBS 

Bucharest/Romania 

16754 

LGIR 

Crete Island/Greece 

37018 

— 

Tuapse/Russia 
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Figure 13. Geographic Locations of the Stations 


The first few lines of sample data are shown in Figure 14. It is a detailed 
data set that includes all necessary parameters to calculate modified refractivity. 


17062 Observations at OOZ 01 Jan 2005 


PRBS 

hPa 

HGHT 

m 

TEMP 

C 

DWPT 

C 

ESLH 

% 

MIXR 

g/kg 

DRCT 

deg 

3KNT 

knot 

THTA 

K 

THTB 

K 

THTV 

K 

1016.0 

39 

8.2 

7.1 

93 

6.26 

40 

12 

280.1 

297.4 

281,1 

1015.0 

47 

8.2 

7.1 

93 

6.26 

80 

7 

280.1 

297.4 

281.2 

1000.0 

167 

7.8 

6.8 

93 

6.23 

80 

13 

280.9 

298.2 

282.0 


Figure 14. Sample Data 
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B. MATLAB CODES 


The MATLAB codes used in this thesis for determining duct statistics are 
based on Stahihut’s (2006) routines that were written to calculate refractivity in 
the Arctic regions. The MATLAB routines were revised to input data from the 
University of Wyoming data sets. The MATLAB code is shown in the Appendix. 

First, all soundings are read by MATLAB. After the modified refractivity is 
calculated for every height, the M gradient is calculated. According to the M 
gradients, the ducting conditions are determined for each profile. The lowest 
three ducts are examined if they exist. According to the duct bottom, surface 
duct existence is checked. Then key duct parameters are recorded for each 
station. Finally, an Excel sheet that includes the duct heights, duct thickness, 
duct gradient and other information related to the station is created. This cycle is 
done for every sounding for a year. All results are collected then saved as an 
Excel (.xls) file. This Excel file is loaded into MATLAB. To find statistics for one 
station the “statistics.m” program (see Appendix, Section D) is run. To find 
statistics for a month the ductformonth.m file is run before “statistics.m.” 

To calculate a modified refractivity at a given height. Equation 3.7 is used. 
Equation 4.1, gives the vapor pressure at a given height where the temperature, 
T, is in Celsius: 


The duct top and duct middle are the local minimum and local maximum of 
M, respectively, found by a switch in the sign of dM/dz. The duct bottom has the 
same value of M as at the duct top. It can be calculated by dropping a straight 
vertical line from duct top to the duct bottom. It can be found by applying simple 
geometry. Equation 4.2 gives the height of duct bottom (Stahihut, 2006). 


z 


= 72 


M^-M* 

M^-M, 


7i) 


(4.2) 
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V. RESULTS 


A. DUCT FREQUENCY STATISTICS 

The objective of this chapter is to evaluate the results gathered from 
MATLAB. In this study it is important to know and understand how ducting 
conditions can effect the performance of communication systems. 

A transmitter or receiver placed in a surface duct will experience 
significant signal enhancement, and long detection range. Omni¬ 
directional sources will experience cylindrical spreading loses, and 
narrow radar beams will view the duct as a waveguide. Existence 
of trapping layer is determined by a temperature increase and/or 
humidity decrease with height. (Cairns-McFeeters, 1992) 

Table 4 shows frequency of occurrence of ducts during the five-year 
period from 2005-2009 for each station. Ducti represents the lowest ducts that 
also include the surface-based ducts in Table 4. Duct2 is the first elevated ducts 
that are higher than the lowest ducts. DuctS is the third highest ducts in M 
profiles. 

The significant result is that the probability of occurrence of a duct is 
relatively high for all locations. If we compare day occurrence with night 
occurrence, it is clear that there is not much differences between them. For 
example, Istanbul’s day occurrence frequency is 50.93 and its night occurrence 
frequency is 50.44. Although this conclusion is also true for other trapping 
layers, some stations have significant occurrence differences between night and 
day occurrences of surface ducts. The surface-based duct occurrence in Ankara, 
Istanbul, and Izmir is higher during day hours than night hours. 
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Table 4. 


Frequency of Ducts 


Present occurrence of ducts during 24 
hours 

Ducti 

Duct2 

Duct3 

Sf. Duct 

17351 Adana 

59.46 

23.38 

6.40 

16.05 

17131 Ankara 

34.22 

5.96 

0.56 

18.77 

17281 Diyarbakir 

37.53 

6.18 

0.56 

21.92 

17240 Isparta 

38.15 

8.29 

0.84 

22.00 

17062 Istanbul 

50.67 

16.54 

4.06 

17.10 

17220 Izmir 

45.59 

11.07 

1.58 

20.16 

16716 Athens 

64.69 

18.94 

3.27 

44.20 

15420 Bucharest 

30.24 

5.74 

0.83 

11.26 

16754 Crete Island 

65.59 

23.81 

5.57 

29.45 

37018 Tuapse 

21.05 

2.33 

0.14 

13.36 


Present occurrence of ducts during day 

Ducti 

Duct2 

Duct3 

Sf. Duct 

17351 Adana 

51.56 

19.34 

4.00 

23.10 

17131 Ankara 

43.53 

9.95 

0.85 

32.00 

17281 Diyarbakir 

34.79 

5.93 

0.81 

22.25 

17240 Isparta 

38.53 

9.73 

1.30 

24.72 

17062 Istanbul 

50.93 

17.45 

4.69 

25.92 

17220 Izmir 

52.85 

14.94 

2.34 

33.35 

16716 Athens 

64.22 

18.75 

2.84 

46.89 

15420 Bucharest 

28.50 

5.13 

0.61 

12.27 

16754 Crete Island 

67.24 

26.55 

6.40 

35.89 

37018 Tuapse 

22.89 

2.53 

0.07 

15.24 


Present occurrence of ducts during night 

Ducti 

Duct2 

Duct3 

Sf. Duct 

17351 Adana 

67.15 

27.32 

8.76 

9.11 

17131 Ankara 

25.10 

2.05 

0.28 

5.82 

17281 Diyarbakir 

40.18 

6.41 

0.33 

21.60 

17240 Isparta 

37.80 

6.90 

0.40 

19.41 

17062 Istanbul 

50.44 

15.65 

3.44 

8.44 

17220 Izmir 

38.46 

7.30 

0.84 

7.30 

16716 Athens 

65.23 

19.19 

3.61 

42.30 

15420 Bucharest 

31.97 

6.35 

1.05 

10.27 

16754 Crete Island 

63.91 

20.87 

4.75 

22.66 

37018 Tuapse 

19.21 

2.13 

0.21 

11.48 
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Table 5. Present Occurrence of First Duct by Month (Surface-based Duct is Included) 



Jan 

Feb 

Mar 

Apr 

May 

Jun 

Jly 

Aug 

Sep 

Oct 

Nov 

Dec 

17351 Adana 

33.44 

29.63 

43.73 

58.87 

68.92 

82.87 

91.52 

94.85 

71.72 

57.89 

42.81 

39.80 

17131 Ankara 

23.53 

25.00 

20.79 

26.44 

43.52 

44.48 

51.97 

51.80 

34.29 

30.00 

28.28 

32.89 

17281 Diyarbakir 

18.72 

27.19 

27.05 

37.45 

56.77 

53.95 

62.81 

55.14 

47.12 

24.57 

23.89 

21.90 

17240 Isparta 

20.20 

16.86 

23.31 

30.98 

45.85 

57.39 

66.54 

67.93 

48.78 

37.46 

25.17 

21.43 

17062 Istanbul 

25.00 

20.73 

29.47 

47.80 

60.98 

72.70 

81.46 

85.99 

59.93 

48.22 

40.94 

33.22 

17220 Izmir 

28.95 

29.71 

32.23 

42.37 

47.04 

61.59 

72.24 

65.25 

56.47 

38.06 

41.75 

31.40 

16716 Athens 

46.49 

44.08 

56.43 

69.51 

72.28 

83.12 

83.06 

83.77 

75.83 

63.30 

57.29 

53.21 

15420 Bucharest 

9.06 

6.16 

14.62 

26.76 

35.20 

51.38 

56.39 

55.19 

37.04 

31.72 

23.08 

14.94 

16754 Crete Island 

44.68 

41.87 

54.31 

62.45 

69.62 

88.70 

84.85 

85.09 

75.90 

66.67 

66.67 

56.79 

37018 Tuapse 

3.42 

8.63 

8.26 

15.86 

27.65 

35.05 

47.37 

40.57 

30.47 

22.69 

12.89 

7.07 


25 




















































































































































J^n 


l-ob Mar Apr May Juti Jly A.ig Sep Ocl Nov Ucc 


■ 17351Arinna 

■ l/lil Ankara 

■ 1/281 Divarbak r 

■ l/240ispafta 

■ 17002 IstHr bul 

■ 177?Olzmi- 


Figure 15. Frequencies of First Duct by Month (Turkey) 
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Figure 16. Frequencies of First Duct by Month (Athens, Bucharest, Crete Island, and Tuapse) 
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Table 5, Figure 15 and Figure 16 all show the distribution of first ducts 
(surface-based ducts are included) by month. It is clear that the ducts are most 
common during summer time. In Adana, there is greater than 90% probability of 
ducts occurring in July and August. The lowest occurrence of first ducts within 
Turkey was in February, at Isparta, with 17 percent of duct occurrence. This 
percentage rate of ducting conditions is significant, indicating that ducting may be 
a factor in all regions of Turkey during all times of the year. The highest 
occurrence of the stations outside of Turkey were Athens and Crete Island, 
where they have almost the same duct frequency trend because they are 
affected by the same climate conditions. 

Comparing these results to AREPS (Advanced Refractive Effects 
Prediction System) climatology clearly shows that the percentage of duct 
occurrences for all stations is higher than AREPS’s. Figure 17 and Figure 18 are 
the sample elevated duct frequency distributions by months for Ankara and 
Isparta from AREPS. 


Ankara /Central, Turkey 



I 1200Z 
I OOOOZ 


“I-1-T 

1 2 3 4 5 

Elevated duct - % occur day and night 


Figure 17. Elevated Duct Percent Occurance for Day and Night in Ankara 
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Isparta, Turkey 



Figure 18. Elevated Duct Percent Occurence for Day and Night in Isparta 

B. DUCT THICKNESS AND HEIGHT STATISTICS 

The following formula, based on experimental results, is used to determine 
the minimum frequency that will be trapped by a duct of given thickness (d, in 
meters): 


/..n= (5.1) 

V m ^ 

f igure 19 represents the equations graphically, while Table 6 summarizes 
duct thickness ranges. Duct thickness statistics of the stations can be evaluated 
according to Table 6. It shows the common frequency ranges and the required 
duct depths to trap a signal in these frequency ranges (Davidson, 2003). 
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Table 6. 


Minimum Duct Thickness Required for Trapping Determined by 
Radar Frequency Band [From (Cairns-McFeeters, 1992)] 


BAND 

Frequency Range 

Duct Thickness (m) 

VHF 

30 MHz-300 MHz 

122-610 

UHF 

300 MHz - 1 GHz 

52-122 

L 

1 GHz-2 GHz 

32-52 

S 

2 GHz-4 GHz 

20-32 

C 

4 GHz-8 GHz 

12-20 

X 

8 GHz - 12 GHz 

9-12 

K 

12 GHz - 40GHz 

4-9 



0 200 60 400120 600 180 800 240 1000 300 

_ Duct Thickness (ft/in) _ 

Figure 19. Minimum Trapping Frequency [From (Davidson, 2003)] 

Table 7 displays the height statistics of ducts and Table 8 shows the 
thickness of the first three and surface ducts. Fleights are above sea level. 
Surface ducts are very shallow for all stations. Note from Table 6 that in many 
cases, only the highest frequency transmissions will be trapped. In addition, 
topographic features such as hills are usually higher than the duct, which also 
would prevent the normally extended ranges associated with ducting situations. 
Although over water and other flat locations microwaves may have extended 
ranges due to the ducts, in most cases these shallow surface ducts are not 
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operationally significant. We can compare duct thicknesses of our locations with 
Table 6 above. There is a reverse relationship between frequency and duct 
thickness. As frequency increases, the required duct thickness to send a signal 
decreases. The mean height above sea level of first ducts in Turkey is around 
8501100 meters for all stations. Bucharest has the highest mean thickness of 
surface ducts at only 74.5 meters. According to Table 6, this would trap radiation 
in the UHF and higher frequency band ranges. If we examine the other layers, 
Adana has the highest mean thickness at 199.9 meters. 
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Table 7. 


Duct Height Statistics 


ADANA (ABOVE SEA LEVEL 25m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

856.2 

1863.6 

2543.6 

59.1 

Lower Quartile (m) 

86.0 

1087.0 

1550.5 

36.0 

Median (m) 

565.0 

1535.5 

2222.5 

45.0 

Higher Quartile (m) 

1346.0 

2504.5 

3552.5 

64.0 

Std. Dev. (m) 

970.1 

1076.6 

1239.3 

34.7 

ANKARA (ABOVE SEA LEVEL 891m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

1824.9 

3169.3 

3266.7 

907.4 

Lower Quartile (m) 

901.0 

2651.0 

2973.0 

900.0 

Median (m) 

920.0 

3186.0 

3381.0 

901.0 

Higher Quartile (m) 

2847.0 

3730.0 

3670.0 

910.0 

Std. Dev. (m) 

1201.3 

851.5 

679.8 

12.8 

DIYARBAKIR (ABOVE SEA LEVEL 675m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

1417.2 

3267.6 

4213.9 

695.8 

Lower Quartile (m) 

685.0 

2362.0 

3589.0 

684.0 

Median (m) 

706.5 

3509.5 

3971.0 

691.0 

Higher Quartile (m) 

1510.0 

4218.0 

4987.5 

700.0 

Std. Dev. (m) 

1293.9 

1413.4 

872.4 

18.4 

ISPARTA (ABOVE SEA LEVEL 997m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

1759.2 

3461.2 

3502.7 

1022.7 

Lower Quartile (m) 

1016.0 

2960.0 

2910.5 

1007.0 

Median (m) 

1036.0 

3507.0 

3285.0 

1016.0 

Higher Quartile (m) 

2706.0 

3949.0 

4030.0 

1027.0 

Std. Dev. (m) 

1122.6 

886.8 

753.7 

19.9 

ISTANBUL (ABOVE SEA LEVEL 39m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

787.9 

1668.0 

2447.5 

63.5 

Lower Quartile (m) 

65.0 

941.5 

1563.0 

48.0 

Median (m) 

419.0 

1485.0 

2373.0 

56.0 

Higher Quartile (m) 

1260.0 

2264.0 

3257.5 

68.0 

Std. Dev. (m) 

927.7 

1003.4 

1094.7 

24.4 

IZMIR (ABOVE SEA LEVEL 29m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

952.9 

2062.5 

2528.6 

51.0 

Lower Quartile (m) 

49.0 

1292.0 

1531.0 

38.0 

Median (m) 

174.0 

1971.0 

2470.0 

46.0 

Higher Quartile (m) 

1718.0 

2718.0 

3404.5 

57.0 

Std. Dev. (m) 

1193.0 

1056.9 

1217.1 

15.5 
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ATHENS (ABOVE SEA LEVEL 15m) 

Duct1 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

505.1 

1864.6 

2747.9 

53.8 

Lower Quartile (m) 

39.0 

743.0 

1928.5 

32.0 

Median (m) 

59.5 

1849.0 

2773.0 

47.0 

Higher Quartiie (m) 

186.0 

2774.0 

3390.0 

65.0 

Std. Dev. (m) 

979.8 

1261.9 

1300.7 

31.3 

BUCHAREST (ABOVE SEA LEVEL 91m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

953.0 

2540.9 

3443.2 

157.6 

Lower Quartiie (m) 

154.0 

1788.0 

2858.0 

139.0 

Median (m) 

253.5 

2474.0 

3407.0 

150.0 

Higher Quartiie (m) 

1585.0 

3161.0 

3953.0 

167.0 

Std. Dev. (m) 

1153.9 

1155.8 

899.7 

33.2 

CRETE ISLAND (ABOVE SEA LEVEL 39m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

761.5 

1515.4 

2428.7 

80.1 

Lower Quartiie (m) 

73.0 

488.0 

1270.0 

56.0 

Median (m) 

189.5 

1193.0 

2126.0 

72.0 

Higher Quartiie (m) 

1238.0 

2142.5 

3311.5 

95.0 

Std. Dev. (m) 

1021.8 

1251.7 

1513.9 

39.0 

TUAPSE (ABOVE SEA LEVEL 95m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Height (m) 

589.2 

1652.5 

1539.5 

133.7 

Lower Quartiie (m) 

116.0 

788.0 

1034.0 

109.5 

Median (m) 

139.0 

1566.0 

1575.0 

121.0 

Higher Quartiie (m) 

782.0 

2207.5 

2045.0 

136.0 

Std. Dev. (m) 

1087.9 

990.8 

608.1 

50.5 
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Table 8. 


Duct Thickness Statistics 


ADANA (ABOVE SEA LEVEL 25m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

147.5 

199.9 

156.6 

58.2 

Lower Quartile (m) 

35.0 

91.0 

78.6 

14.0 

Median (m) 

104.3 

155.5 

119.3 

26.0 

Higher Quartile (m) 

212.2 

277.6 

184.9 

139.0 

Std. Dev. (m) 

140.4 

144.8 

121.2 

74.9 

ANKARA (ABOVE SEA LEVEL 891m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

52.2 

95.1 

69.7 

16.4 

Lower Quartile (m) 

10.0 

59.0 

41.1 

9.0 

Median (m) 

28.0 

87.9 

62.2 

10.0 

Higher Quartile (m) 

87.3 

124.7 

88.4 

58.0 

Std. Dev. (m) 

52.5 

49.7 

34.3 

12.8 

DIYARBAKIR (ABOVE SEA LEVEL 675m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

45.2 

81.8 

79.3 

22.6 

Lower Quartile (m) 

10.0 

46.3 

47.8 

9.0 

Median (m) 

26.0 

75.0 

80.6 

17.0 

Higher Quartile (m) 

68.8 

113.9 

112.8 

52.1 

Std. Dev. (m) 

45.0 

48.6 

36.2 

22.0 

ISPARTA (ABOVE SEA LEVEL 997m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

51.8 

90.3 

66.3 

28.4 

Lower Quartile (m) 

19.0 

60.7 

39.4 

10.0 

Median (m) 

36.2 

83.4 

50.6 

19.0 

Higher Quartile (m) 

75.9 

116.5 

92.5 

57.0 

Std. Dev. (m) 

45.7 

41.5 

40.3 

24.5 

ISTANBUL (ABOVE SEA LEVEL 39m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

116.4 

153.7 

122.4 

41.8 

Lower Quartile (m) 

24.0 

85.9 

60.3 

9.0 

Median (m) 

84.0 

133.5 

104.5 

18.0 

Higher Quartile (m) 

170.0 

191.1 

166.8 

112.6 

Std. Dev. (m) 

113.3 

100.2 

79.7 

69.5 

IZMIR (ABOVE SEA LEVEL 29m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

70.4 

121.7 

118.0 

26.6 

Lower Quartile (m) 

17.0 

76.9 

53.9 

9.0 

Median (m) 

44.0 

112.6 

90.7 

17.0 

Higher Quartile (m) 

105.6 

158.5 

160.7 

71.9 

Std. Dev. (m) 

67.3 

66.6 

80.9 

31.9 
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ATHENS (ABOVE SEA LEVEL 15m) 

Duct1 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

59.3 

116.9 

99.3 

42.0 

Lower Quartile (m) 

24.0 

69.4 

63.0 

17.0 

Median (m) 

40.0 

106.1 

86.9 

32.0 

Higher Quartile (m) 

81.6 

154.1 

123.6 

68.9 

Std. Dev. (m) 

54.5 

67.8 

53.8 

37.4 

BUCHAREST (ABOVE SEA LEVEL 91m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

93.8 

126.5 

87.7 

74.5 

Lower Quartile (m) 

53.0 

79.6 

45.6 

47.2 

Median (m) 

73.8 

111.1 

77.7 

61.0 

Higher Quartile (m) 

121.0 

167.0 

131.0 

109.0 

Std. Dev. (m) 

60.6 

67.0 

52.7 

43.5 

CRETE ISLAND (ABOVE SEA LEVEL 675m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

96.5 

138.5 

107.0 

51.8 

Lower Quartile (m) 

28.6 

71.0 

50.3 

19.0 

Median (m) 

65.0 

124.8 

93.7 

34.0 

Higher Quartile (m) 

139.7 

188.3 

139.5 

106.0 

Std. Dev. (m) 

90.2 

87.6 

82.6 

59.0 

TUAPSE (ABOVE SEA LEVEL 997m) 

Ducti 

Duct2 

Duct3 

Sf. Duct 

Mean Thickness (m) 

76.2 

135.8 

110.9 

39.8 

Lower Quartile (m) 

19.0 

84.0 

55.8 

15.0 

Median (m) 

38.0 

106.5 

105.9 

26.0 

Higher Quartile (m) 

96.9 

156.2 

165.9 

68.0 

Std. Dev. (m) 

94.0 

86.8 

75.6 

51.7 
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c. 


CORRELATIONS WITH METEOROLOGICAL PARAMETERS 


To provide guidance on what factors may be related to duct occurrence, 
surface pressure, wind speed, and wind direction were plotted against duct 
height to determine if there are any correlations that could be used to make more 
accurate predictions of duct occurrence. Although these meteorological 
parameters are not directly related to ducting (which depends primarily on vertical 
variations in humidity), the hypothesis is that these parameters may affect the 
humidity profiles which in turn would be related to ducting probabilities. For 
example, wind speed may affect atmospheric mixing of water vapor (humidity) 
which results in changes in ducting statistics. Wind direction could be correlated 
with moisture advection (horizontal movement), especially in coastal regions 
where it would be expected than air coming from over the water would have 
different humidity, and hence ducting characteristics. Higher surface pressure is 
associated with sinking air motion, which is often associated with the formation of 
an elevated temperature inversion (increase in temperature with height), and a 
strong decrease in humidity going up through the inversion. It is well known that 
in most locations ducting is correlated with high pressure. 

Scatter plots of these atmospheric parameters and duct height were 
created for all stations. Only those stations with significant correlations are 
displayed in this section. The most significant conclusion for wind speed and 
wind direction is that they do not appear to have a significant effect on duct 
heights for most stations, though some exceptions are shown below. It can also 
be seen that, from plots for some stations, there is a positive relation between 
duct height and surface pressure. This is somewhat unexpected because the 
sinking air associated with high pressure was expected to push the duct heights 
lower. 
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Mtop-Surface Pressure, (flhere is dud, Adana, Above Sea Level: 25m 



Figure 20. Height (m) and Surface Pressure (milibar) Correlation for Adana (if 

a Duct Presents) 

Figure 20 shows that a higher pressure was associated with lower ducts at 
Adana. Isparta (Figure 22) also has the same pattern. If we examine Figure 21, 
which is for surface ducts only, it is clearly determined that if the surface pressure 
increases, the duct height also increases. This may be because the shallower 
ducts are associated with shallow temperature inversions caused by radiational 
cooling (which does not require high pressure) while the deeper ducts may be 
more associated with the elevated inversions caused by sinking air. 
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Mtop-Surface Pressure, if there is surface duct, Adana, Above Sea Level: 25m 



Figure 21. Height (m) and Surface Pressure (milibar) Correlation for Adana (if 

only a Surface Duct Presents) 
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Figure 22. Height (m) and Surface Pressure (milibar) Correlation for Isparta (if 

a Duct Presents) 
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Wop-Wind Speed, if there is duct. Isparta, Above Sea Level; 997m 



Figure 23. Height (m) and Wind Speed (knot) Correlation for Isparta (if a Duct 

Presents) 

For Isparta, there is a positive relationship between wind speed and duct 
height if there is any type of duct (Figure 23). 

Figure 24 and Figure 28 are other examples of a positive relationship 
between surface pressure and duct height, if any duct is present or if only a 
surface duct exists. Figure 26 shows that if the wind direction switches form east 
through north in Istanbul, it is possible to see more numbers of ducts. In 
Bucharest, if the wind speed increases the duct height also increases (Figure 
27). 
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Miop-Sutface Pressure. IfIhere is duel, Istanbul, Above Sea Level 39m 



Figure 24. Height (m) and Surface Pressure (milibar) Correlation for Istanbul (if 

a Duct Presents) 


Mtop-Surface Pressure, if there is surface duct. Istanbul. Above Sea Level 39m 



Figure 25. Height (m) and Surface Pressure (milibar) Correlation for Istanbul (if 

Only a Surface Duct Presents) 
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Figure 26. Height (m) and Wind Direction (degree) Correlation for Istanbul (if a 

Duct Presents) 

fvltop-Wind Speed, if there is dud, Bucharest, Above Sea Level: 91m 



Wind Speed 


Figure 27. Height (m) and Wind Speed (knot) Correlation for Bucharest (if a 

Duct Presents) 
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Ml Op-Surface Pressure, If there is duct. Crete Island, Above Sea Level 39m 



Figure 28. Height (m) and Surface Pressure (milibar) Correlation for Crete 

Island (if a Duct Presents) 
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VI. CONCLUSION 


The purpose of this thesis was to characterize ducting climatology within 
Turkey and its surrounding areas. Almost 36,500 lines of data, representing 
different elevations, stations and times, were processed. Refractivity profile 
statistics for each station were calculated. The general trend of duct occurrences 
and the relationship between duct height and other parameters were quantified. 

According to the results, the probability of duct occurrence is higher in 
summer months. Wind speed and wind direction do not have much effect on 
ducting for most stations. There is a positive relationship between surface 
pressure and duct height at some stations. This condition was not an expected 
situation. A likely explanation is that very shallow ducts, associated with surface 
radiational cooling or other local effects, are often present during both high and 
low pressure situations, but the higher ducts are associated more with higher 
pressure. Therefore, when considering all ducts, the average height is higher for 
high-pressure situations. Although there are many ducts, they are not so 
important for VHP and lower frequencies because their thickness is not enough 
to send a signal in VHP and lower frequency ranges. These ducts can be 
important for systems that use a microwave frequency range (UPlP and EPIS). 
Por example, radar uses microwave radiation to detect the range of remote 
objects. 

The plots, tables, and figures provided in this thesis can be a guide for 
Electronic Warfare Systems and radio operators and their commanders. By using 
knowledge of the environment and its affect on radio transmissions, we can gain 
an advantage over adversaries who do not consider this information in planning 
and carrying out military operations. It can also be helpful for civilian use of radio 
transmissions. The most significant threat is losing communication links. 
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Knowing the technical information about an EW system and having the ability to 
predict atmospheric conditions can greatly enhance the communication 
capabilities of Turkey and its allies. 
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APPENDIX 


A. MATLAB PROGRAM TO READ DATA 

function loadsnd{filename) 
fid_in=fopen{filename, 'r' ) ; 
if fid_in > 0 

sprintf {' Sounding file opened: %s filename) 

end 

tic 

% counter for all data 
bigcounter = 0; 

%This loop for all year 
while 1 

% Skip first five rows header 
for ihead = 1:5; 

line = fgetl{fid_in); 

end 

"5 

%initialize a line counter 
iline=0; 

'o 

%This loop reads each line 
while 1 

%Read each line as a character string 
line=fgetl{fid_in) ; 

%This loop finds missing value for dwpt 
while line (22:28) == ' ' 

line=fgetl{fid_in); 

end 

%This is statement for reading headers 
if line (1:7) == 'Station' 
line=fgetl{fid_in); 
line=fgetl{fid_in); 
station = str2num{line{46:50)); 
line=fgetl{fid_in); 
date=str2num{line{4 6:51)); 
time=str2num(line(53:5 6)) ; 
line=fgetl{fid_in); 

%This loop skip other lines 
while ~isempty(line) %line(44) == ':' 
line=fgetl{fid_in); 

%Check for the effective "end-of-file" 
if line == -1 
break 

end 

end 

end 

%Check for the end of one data set 
if isempty(line) 
break 

end 
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%Check: for the effective "end-of-file" 
if line == -1 
break 

end 

%Increment line counter. This is the index for each data array 
iline = iline+1; 

%Reads data from file 
pres{iline) = str2num{line{1:7)) ; 
disp(pres(iline)); 

hght(iline) = str2num(line(8:14)); 
tmpc(iline) = str2num(line(15:21)) ; 
dwpc(iline) = str2num(line(22:28)); 
relh(iline) = str2num(line(29:35)); 
mixr(iline) = str2num(line(36:42)); 


if line(48:56) 
drct(iline) 
sknt(iline) 

else 

drct(iline) 
sknt(iline) 

end 


NaN; 

NaN; 

str2num(line(4 7:4 9)); 
str2num(line(54:56)); 


end % This is the end of the WHILE Loop 


% Calculate theta & q 

[thta,q] = theta_q(pres, tmpC, relh) ; 

% Calculate M and N 

[M,N] = m_n_profile(pres,tmpC, relh, hght) ; 

% Calculates the Modified Refractivity gradient 

Mgrad=NaN; 

for n=l:iline-l ; 

Mgrad(n,1)=(M(n+1) - M(n)) /(hght(n+l) - hght(n)); 

end 

bigcounter = bigcounter + 1; 

% This portion finds the first three ducts, if 

% they exist 

ducttop=[0 0 0]; 

ductmiddle=[0 0 0]; 

ductbottom=[0 0 0]; 

Ttop=NaN; 

Dewtop=NaN; 

Presstop=NaN; 

Mtop=NaN; 

Tmid=NaN; 

Dewmid=NaN; 

Pressmid=NaN; 

Mmid=NaN; 

Mtop2=NaN; 

Mmid2=NaN; 

Mtop3=NaN; 

Mmid3=NaN; 
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% Find tops and mids of ducts.... 

topcounter=0; 

midcounter=0 ; 

topindex=[0 0 0]; 

midindex=[0 0 0]; 

for m=2:length(Mgrad) 

if {(Mgrad(m-l) < 0) && (Mgrad 

topcounter=topcounter+l ; 
topindex(topcounter)=m; 


end 

if ((Mgrad(m-l) > 0) && (Mgrad 

midcounter=midcounter+l; 
midindex(midcounter)=m; 


end 


(m) 


(m) 


> 0 ) ) 


< 0 ) ) 


end 

% Align these depending on case. 

% If the gradient starts negative from the surface, 

% the subsequent indexing is different that for the 
% case if the gradient starts positive from the surface, 
if (Mgrad(l) < 0) 

if (topindex(l) > 0) 

ducttop(1)=hght(topindex(1)) ; 
ductmiddle(1)=hght(1) ; 
ductbottom(1)=hght(1); 

end 

if (topindex(2) > 0) && (midindex(1) > 0) 

ducttop(2)=hght(topindex(2)) ; 
ductmiddle(2)=hght(midindex(1)) ; 
if M(midindex(1)-1) < M(topindex(2)) 

ductbottom(2)=ductmiddle(2) - ... 

(M(midindex(1))-M(topindex(2 )))/.. . 

(M(midindex(1))-M(midindex(1) -1) ) * . . . 

(hght(midindex(1)) - hght(midindex(1)-1)); 

elseif M(midindex(1)-1) > M(topindex(2))&& M(midindex(1) 

2) >M(midindex(1)-1) 

ductbottom(2)=hght(midindex(1)-1); 
elseif M(midindex(1)-2) > M(topindex(2))&& M(midindex(1) 

3) >M(midindex(1)-2) 

ductbottom(2)=hght(midindex(1)-2); 
elseif M(midindex(1)-1) > M(topindex(2)) 

counterindex = 1; 

while M(midindex(1)-counterindex)>M(topindex(2))& & . . 
(midindex(1)-counterindex)~=1 
counterindex = counterindex+1; 

end 

ductbottom(2)=hght(midindex(1)-counterindex+2 ) - ... 

(M(midindex(1)-counterindex+2)-M(topindex(2)))/ . 
(M(midindex(1)-counterindex+2)-M(midindex(1)- 

counterindex )) * ... 

(hght(midindex(1)-counterindex+2) - 

hght(midindex(1)-counterindex) ) ; 
end 

Mtop2=M(topindex(2)) ; 

Mmid2=M(midindex(1)) ; 

end 


47 



if (topindex{3) > 0) && {midindex{2) > 0) 

ducttop(3)=hght(topindex(3)); 
ductmiddle(3)=hght(midindex{2)) ; 
if M{midindex{2)-1) < M(topindex{3)) 
ductbottom{3)=ductmiddle{3) - ... 

(M{midindex{2))-M(topindex{3 )))/.. . 

(M{midindex{2))-M{midindex{2)-1) ) * . . . 

(hght{midindex{2)) - hght{midindex{2)-1)); 

elseif M{midindex{2)-1) > M(topindex{3)) 
counterindex = 1; 

while M{midindex{2)-counterindex)>M(topindex{3) ) 
counterindex = counterindex+1; 

end 

ductbottom{3)=hght{midindex{2)-counterindex+2 ) - .. 

(M{midindex{2)-counterindex+2)-M(topindex{3)))/ 
(M{midindex{2)-counterindex+2)-M{midindex(2)- 

counterindex))* ... 

(hght{midindex{2)-counterindex+2) - 

hght{midindex{2)-counterindex)) ; 
end 

Mtop3=M(topindex{3)) ; 

Mmid3=M{midindex{2)) ; 

end 

Ttop=tmpc(topindex{1)) ; 

Dewtop=dwpc(topindex{1)) ; 

Presstop=pres(topindex{1)); 

Mtop=M(topindex{1)) ; 

Tmid=tmpc{1) ; 

Dewmid=dwpc{1); 

Pressmid=pres(1); 

Mmid=M{l); 

end 

% These are for the cases when the gradient starts 
% positive, and there is a duct present somewhere in 
% the profile 

if {Mgrad{l) > 0) & {topindex{l) > 0) 

if {topindex{l) > 0) & {midindex{l) >0) & ... 

(M(topindex{1)) < M{1)) 

ducttop{1)=hght(topindex{1)) ; 
ductmiddle{1)=hght{midindex{1)) ; 
ductbottom{1)=hght{1) ; 

end 

if {topindex{l) > 0) & {midindex{l) >0) & ... 

(M(topindex{1)) > M{1)) 

ducttop{1)=hght(topindex{1)) ; 
ductmiddle{1)=hght{midindex{1) ) ; 
if M{midindex{1)-1)<M(topindex{1)) 

ductbottom{1)=ductmiddle{1) - ... 

(M{midindex{1))-M(topindex(1 )))/.. . 

(M{midindex{1))-M{midindex{1)-1))* ... 

(hght{midindex{1)) - hght{midindex{1)-1)); 

elseif M{midindex{1)-1)>M(topindex{1)) 
counterindex = 1; 
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while M{midindex(1)-counterindex)>M(topindex(1)) 
counterindex = counterindex+1; 

end 

ductbottom(1)=hght(midindex(1)-counterindex+2 ) - .. 

(M(midindex{1)-counterindex+2)-M(topindex {!)))/ 
(M{midindex{1)-counterindex+2)-M(midindex(1)- 

counterindex))* ... 

{(hght(midindex(1)-counterindex+2)) - 

hght{midindex{1)-counterindex) ) ; 
end 

end 

if (topindex(2) > 0) & {midindex(2) > 0) 

ducttop{2)=hght(topindex{2)); 
ductmiddle{2)=hght{midindex(2 ) ) ; 
ductbottom{2)=ductmiddle{2) - ... 

(M(midindex(2))-M(topindex(2 )))/.. . 

(M{midindex{2))-M(midindex{2)-1) ) * . . . 

(hght{midindex{2)) - hght{midindex(2)-1)); 

if M{midindex{2)-1)<M(topindex{2)) 

ductbottom{2)=ductmiddle{2) - ... 

(M(midindex(2))-M(topindex{2 )))/.. . 

(M{midindex{2))-M(midindex(2)-1))* ... 

(hght{midindex{2)) - hght{midindex{2)-1)); 

elseif M{midindex{2)-1)>M(topindex{2)) 
counterindex = 1; 

while M(midindex{2)-counterindex)>M(topindex{2))& & . 
(midindex(2)-counterindex) ~=1 
counterindex = counterindex+1; 

end 

ductbottom(2)=hght{midindex(2)-counterindex+2 ) - .. 

(M(midindex(2)-counterindex+2)-M(topindex(2)))/ 
(M{midindex{2)-counterindex+2)-M{midindex(2)- 

counterindex))* ... 

{(hght(midindex(2)-counterindex+2)) - 

hght{midindex{2)-counterindex)) ; 
end 

Mtop2=M(topindex{2)) ; 

Mmid2=M{midindex{2)) ; 

end 

if (topindex(3) > 0) & (midindex(3) > 0) 
ducttop(3)=hght(topindex(3)); 
ductmiddle(3)=hght{midindex(3)) ; 
ductbottom{3)=ductmiddle{3) - ... 

(M{midindex{3))-M(topindex{3 )))/.. . 

(M(midindex(3))-M{midindex{3)-1) ) * . . . 

(hght(midindex(3)) - hght{midindex(3)-1)); 
if M{midindex(3)-1)<M(topindex(3)) 

ductbottom{3)=ductmiddle{3) - ... 

(M{midindex{3))-M(topindex{3 )))/.. . 

(M(midindex(3))-M{midindex{3)-1))* ... 

(hght(midindex{3)) - hght(midindex{3)-1)); 
elseif M(midindex(3)-1)>M(topindex(3)) 
counterindex = 1; 
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while M{midindex(3)-counterindex)>M(topindex(3)) 
counterindex = counterindex+1; 

end 

ductbottom(3)=hght(midindex(3)-counterindex+2 ) - .. 

(M(midindex{3)-counterindex+2)-M(topindex(3)))/ 
(M{midindex(3)-counterindex+2)-M{midindex(3)- 

counterindex))* ... 

{(hght{midindex{3)-counterindex+2))- 
hght{midindex(3)-counterindex) ) ; 
end 

Mtop3=M(topindex{3)) ; 

Mmid3=M{midindex{3)) ; 

end 

Ttop=tmpc(topindex{1)) ; 

Dewtop=dwpc(topindex{1)) ; 

Presstop=pres(topindex{1)); 

Mtop=M(topindex(1)) ; 

Tmid=tmpc{midindex{1)) ; 

Dewmid=dwpc(midindex{1) ) ; 

Pressmid=pres{midindex(1) ) ; 

Mmid=M{midindex(1)) ; 

if topindex(2)~=0 && midindex{2)~=0 
Mtop2=M(topindex{2)) ; 

Mmid2=M{midindex{2)) ; 

else 

Mtop2=NaN; 

Mmid2=NaN; 

end 

if topindex{3)~=0 && midindex{3)~=0 
Mtop3=M(topindex{3)); 

Mmid3=M(midindex{3)) ; 

else 

Mtop3=NaN; 

Mmid3=NaN; 

end 

end 

if ductbottom{1) < 0 

ductbottom{1)=ductmiddle(1)-[ducttop{1)-ductmiddle(1)] ; 

end 

if ductbottom(2) < 0 

ductbottom{2)=ductmiddle{2)-[ducttop{2)-ductmiddle(2)] ; 

end 

if ductbottom(1) < 0 

ductbottom{1) =0; 

end 

if ductbottom{2) < 0 

ductbottom{2) =0; 

end 
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% This will assign a true/false value if 
% the lowest duct is "attached" to or 
% "elevated" from the surface 

if [ducttop(l) > 0] & [ductbottom{1) > hght(l)] 
sfcduct=0; 

elseif {[ducttop{l) > 0] & [ductbottom{1) == hght{l)]) |... 
{[ducttop{l) >0] & [ductbottom{1) == 0]) 
sfcduct=l; 

else 

sfcduct=NaN; 

end 

% Assigns needed data to the Plotdata matrix 
Plotdata(bigcounter,l)=station; 

Plotdata(bigcounter,2)=date; 

Plotdata(bigcounter,3)=time; 

Plotdata(bigcounter,4)=Ttop; 

Plotdata(bigcounter,5)=Dewtop; 

Plotdata(bigcounter,6)=ducttop(1); 

Plotdata(bigcounter,7)=ductmiddle(1) ; 

Plotdata(bigcounter,8)=ductbottom(1); 

Plotdata(bigcounter,9)=ducttop(2) ; 

Plotdata(bigcounter,10)=ductmiddle(2); 

Plotdata(bigcounter,11)=ductbottom(2); 

Plotdata(bigcounter,12)=ducttop(3); 

Plotdata(bigcounter,13)=ductmiddle(3); 

Plotdata(bigcounter,14)=ductbottom(3); 

Plotdata(bigcounter,15)=sfcduct; 

Plotdata(bigcounter,16)=Presstop; 

Plotdata(bigcounter,17)=Pressmid; 

Plotdata(bigcounter,18)=Tmid; 

Plotdata(bigcounter,19)=Dewmid; 

Plotdata(bigcounter,20)=Mtop; 

Plotdata(bigcounter,21)=Mmid; 

Plotdata(bigcounter,22)=pres(1) ; 

Plotdata(bigcounter,23)=drct(1) ; 

Plotdata(bigcounter,24)=sknt(1) ; 

Plotdata(bigcounter,25)=ducttop(1)-ductbottom(1); 

Plotdata(bigcounter,26)=ducttop(2)-ductbottom(2); 

Plotdata(bigcounter, 21 ) =ducttop(3)-ductbottom(3); 

Plotdata(bigcounter,2 8)=Mtop2; 

Plotdata(bigcounter,2 9)=Mmid2; 

Plotdata(bigcounter,30)=Mtop3; 

Plotdata(bigcounter,31)=Mmid3; 

Plotdata(bigcounter,32)=Mmid-Mtop; 

Plotdata(bigcounter,33)=Mmid2-Mtop2; 

Plotdata(bigcounter,34)=Mmid3-Mtop3; 

Plotdata(bigcounter,35)=(Mmid-Mtop)/(ducttop(1)-ductmiddle(1)); 
Plotdata(bigcounter,36) = (Mmid2-Mtop2)/(ducttop(2)-ductmiddle(2)) ; 
Plotdata(bigcounter,37) = (Mmid3-Mtop3)/(ducttop(3)-ductmiddle(3)) ; 
Plotdata(bigcounter,38)=M(1); 

if Plotdata(bigcounter,15)==1 %finds M deficit 
Plotdata(bigcounter,39) = M(l)-Mtop; 

else 

Plotdata(bigcounter,39) = NaN; 

end 
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end 


%Check for the effective "end-of-file" 
if line == -1 
break 

end 


fclose{fid_in); % - end of soundrng "decoder" portron 


% prepare the name for the ".mat" fi 
% Logic: find the "/"s and so 

name. 

"o 

indexl=find(filename == 
if length{indexl) == 0 
first=l; 

else 

iend=length{indexl) ; 
first=indexl(lend) +1; 

name 

end 

index2 = find(filename == ' . ' ) ; 
last=index2-l ; 

matfile=filename(1,first:last) ; 

%clear indexl index2 first last lend 

"o 

eval{[ 'save ' matfile] ) 

"o 

disp{'".mat" file has been written ti 
toe 


e 

e can use the first part of the 

% find the "/" in the name 
% if no "/" in name 

% the last "/" 

% first character in matfile 

% use characters before the "." 
ans filename jo iline 

your current directory.') 
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B. MATLAB PROGRAM TO CALCULATE MODIFIED REFRACTIVITY 


function [M,N] = m_n_profile(pres,tmpc, relh, hght) 

g, 

o 

% Purpose: Calculates the vertical profiles of Refractive Index, N, and 
% the Modified Refractive Index, M, for input vertical 

% profiles of Pressure, Temperature Relative Humidity and 


Height. 


"O 

% Input: 

pres 

"O 

tmpC 

"o 

relh 

"o 

o. 

hght 

o 

% Output: 
units) 

M = 1 


N = Refractive Index (dimensionless refractivity units) 
Local Variables: 


% tmpk = air temperature (Kelvin) 

% e_s = saturation vapor pressure (millibars) 

% ee = vapor pressure (millibars) 

"6 

% References 

% Bean and Dutton, 1968, Equation 1.16. ... for N 

% Patterson, et al. (1994), Equation 5, p. 9. ... for M 

% Bolton, Monthly Weather Review, 1980. ... for saturation vapor 

pres 

% Huschke, Glossary of Meteorology, 1959, p. 477. ... for vapor 

pressure 

"o 

% History 

% Version 1.0 17 July 2001 

% Mary S. Jordan 

% Dept, of Meteorology, Naval Postgraduate School, Monterey, CA 


% ... convert to temperature to Kelvin 

tmpk = tmpc+273.15; 

% ... compute saturation vapor pressure, e_s 

e_s = 6.112*exp((17.67 .* tmpc) ./(tmpc+243.5)); 

% ... compute vapor pressure, ee 

ee = (relh ./ 100.) .* e_s; 

% ... compute Refractive Index, N, and Modified Refractive Index, M 

"o 

N = (77.6*pres ./tmpk) - (5.6*ee ./tmpk) + (3.75e5*ee . / (tmpk . "'2 ) ) ; 

M = N + 157*hght./1000; 


end of function 


53 





C. MATLAB POGRAM TO CALCULATE POTENTIAL TEMPERATURE 
AND SPECIFIC HUMIDITY 


function [thta,q] = theta_q(pres,tmpc, relh) 

"o 

% function [thta,q] = theta_q(pres,tmpc,relh) 

"O 

% written by: Mary Jordan, NPS Meteorology Dept, 10/9/96 

"O 

% Purpose: Calculates Potential Temperature (thta) and 
% Specific Humidity (q) 

"O 

% Reference: Atmospheric Science by Wallace & Hobbs, 1977, Academic 
Press 


% Input: pres = Pressure (mb) 

% tmpC = Air Temperature (C) 

% relh = Relative Humidity (%) 

"O 

% Output: thta = Potential Temperature (K) 
% q = Specific Humidity (g/kg) 

"6 

"o 


"o 

% define local constants for thermodynamic equations: 
p0=1000.;R=287.;cp=1004.;L=2.5e6; 

"6 

tmpk = tmpc+273.155; % convert to tmpc to Kelvin 

"o 

e_s = 6.1078*exp((17.26939.*tmpc)./(tmpc+237.3)); % saturation vapor 

pressure 

ee = (relh ./ 100.) .* e_s; % vapor pressure 


w = 0.622 * (ee ./ (pres - ee)); % mixing ratio (kg/kg) 

q_kg = w ./ (1.0 + w); % specific humidity 

(kg/kg) 

"o 

q=q_kg*1000; % specific humidity (g/kg) 

"6 

thta = tmpk .* ((pO ./ pres).^(R/cp)); % Potential Temperature 

(K) 


END 
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D. 


MATLAB CODE TO FIND STATISTICS 


x=length(data{:,1)); % number of rows 

templ=0; %temperorY variables 

temp2=0; 

temp3=0 ; 

tmpl=0; 

tmp2=0; 

tmp3=0; 

tempminl=0 ; 

tempmin2=0; 

tempmin3=0 ; 

temp=0; 

%finds elevation and sf duct frequencies 
counter = 1; 
ductsl=0; 
for i=l: X 

if data(counter,6)>0 
ductsl=ductsl+l; 

end 

counter=counter+l ; 

end 

percentlduct=(ductsl/x) *100; 
table=[]; 

table(1,1)=percentlduct; 
counter = 1; 
ducts2=0; 

for i=l:X 

if data(counter, 9)>0 
ducts2=ducts2 + l ; 

end 

counter=counter+1 ; 

end 

percent2duct=(ducts2/x) *100; 
table(1,2)=percent2duct ; 
counter = 1; 
ducts3=0; 

for i=l:X 

if data(counter, 12) >0 
ducts3=ducts3 + l ; 

end 

counter=counter+l ; 

end 

percent3duct=(ducts3/x) *100; 
table(1,3)=percent3duct ; 
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%finds frequency of day ductl 

counter=l; 

dayductl=0; 

day=0; 

for i=l: X 

if data(counter,3)==1200 
day=day+l; 

end 

if data(counter,3)==1200 && data(counter,6)>0 
dayductl=dayductl+l; 

end 

counter=counter+l ; 

end 

percentlductday=(dayductl/day)*100; 
table(2,1)=percentlductday; 

%finds frequencu of night ductl 

counter=l; 

nightductl=0; 

night=0; 

for i=l: X 

if data(counter,3)==0 
night=night+l; 

end 

if data(counter,3)==0 && data(counter,6)>0 
nightductl=nightduct1+1; 

end 

counter=counter+l ; 

end 

percentlductnight=(nightductl/night)*100; 
table(3,1)=percentlductnight; 

%finds frequency of day duct2 

counter=l; 

dayduct2=0; 

day=0; 

for i=l: X 

if data(counter,3)==1200 
day=day+l; 

end 

if data(counter,3)==1200 && data(counter,9)>0 
dayduct2=dayduct2+l; 

end 

counter=counter+l ; 

end 

percent2ductday=(dayduct2/day)*100; 
table(2,2)=percent2ductday; 

%finds frequency of night duct2 

counter=l; 

nightduct2=0; 

night=0; 

for i=l: X 

if data(counter,3)==0 
night=night+l; 
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end 

if data(counter,3)==0 && data(counter,9)>0 
nightduct2=nightduct2+l; 

end 

counter=counter+l ; 

end 

percent2ductnight={nightduct2/night)*100; 
table{3,2)=percent2ductnight; 

%finds frequency of day duct3 

counter=l; 

dayduct3=0 ; 

day=0; 

for i=l: X 

if data(counter,3)==1200 
day=day+l; 

end 

if data(counter,3)==1200 && data(counter,12)>0 
dayduct3=dayduct3+l; 

end 

counter=counter+1 ; 

end 

percent3ductday=(dayduct3/day)*100; 
table(2,3)=percent3ductday; 

%finds frequency of night duct3 

counter=l; 

nightduct3=0; 

night=0; 

for 1=1: X 

if data(counter,3)==0 
night=night+l; 

end 

if data(counter,3)==0 && data(counter,12)>0 
nightduct3=nightduct3+l; 

end 

counter=counter+1 ; 

end 

percent3ductnight=(nightduct3/night)*100; 
table(3,3)=percent3ductnight; 

%find mean height of first duct 
counter = 1; 

j=i; 

for 1=1:X 

if data(counter, 6) >0 

tempi(j)=data(counter, 6) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(4,1)=mean(tempi) ; 
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%find mean height of second duct 
counter = 1; 

j=i; 

for i=l:X 

if data(counter, 9)>0 

temp2{j)=data{counter, 9) ; 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(4,2)=mean{temp2) ; 

%find mean height of third duct 
counter = 1; 

j=i; 

for i=l:X 

if data{counter, 12) >0 

temp3{j)=data{counter, 12) ; 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(4,3)=mean{temp3); 

% find medians of ducts 
table(6,1)=median(tempi) ; 
table(6,2)=median(temp2) ; 
table(6,3)=median(temp3) ; 

% compute 25th percentile (first quartile) for ductl 
table(5,1) = median(tempi(find(templ<median(tempi)))) 

% compute 75th percentile (third quartile) for ductl 
table(7,1) = median(tempi(find(templ>median(tempi)))) 

% compute 25th percentile (first quartile) for duct2 
table(5,2) = median(temp2(find(temp2<median(temp2)))) 

% compute 75th percentile (third quartile) for duct2 
table(7,2) = median(temp2(find(temp2>median(temp2)))) 

% compute 25th percentile (first quartile) for ductl 
table(5,3) = median(temp3(find(temp3<median(temp3)))) 

% compute 75th percentile (third quartile) for ductl 
table(7,3) = median(temp3(find(temp3>median(temp3)))) 

% standard deviations of ducts 

table(8,l)=std(tempi); 
table(8,2)=std(temp2) ; 
table(8,3)=std(temp3) ; 


% find mean, median and std. dev. of thickness 
% ductl 
counter = 1; 

j=i; 

templ = 0 ; 
for i=l:X 

if data(counter, 25) >0 

tempi{j)=data{counter,25) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(9,1)=mean(tempi) ; 
table (11,1)=median(tempi) ; 
table(13,l)=std(tempi) ; 

% compute 25th percentile (first quartile) for ductl 
table(10,1) = median(tempi(find(templ<median(tempi)))) ; 

% compute 75th percentile (third quartile) for ductl 
table(12,1) = median(tempi(find(templ>median(tempi)))); 

%duct2 

counter = 1; 

j=i; 

temp2 = 0 ; 
for i=l:X 

if data(counter, 26) >0 

temp2(j)=data(counter, 2 6); 

j=j+i; 

end 

counter=counter+l ; 

end 

table(9,2)=mean(temp2) ; 
table (11,2)=median(temp2) ; 
table(13,2)=std(temp2) ; 

% compute 25th percentile (first quartile) for duct2 
table(10,2) = median(temp2(find(temp2<median(temp2)))); 

% compute 75th percentile (third quartile) for duct2 
table(12,2) = median(temp2(find(temp2>median(temp2)))) ; 

%duct3 

counter = 1; 

j=i; 

temp3=0 ; 
for i=l:X 

if data(counter, 27) >0 

temp3(j)=data(counter, 27) ; 

j=j+i; 

end 

counter=counter+l ; 

end 
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table(9,3)=mean{temp3); 
table (11,3)=median(temp3); 
table{13,3)=std{temp3); 

% compute 25th percentile (first quartile) for duct3 
table(10,3) = median{temp3{find{temp3<median{temp3)))) ; 

% compute 75th percentile (third quartile) for duct3 
table(12,3) = median(temp3(find(temp3>median(temp3)))); 

% frequency of surface duct 
counter=l; 

j=i; 

temp=0 ; 
for i=l:X 

if data(counter,7)==data(counter,8) && data(counter,7)>0 

temp(j)=data(counter, 6) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

a=length(temp); 

table(1,4)=(a/x)*100; %sf duct in all day 
table(4,4)=mean(temp) ; 
table (6,4)=median(temp) ; 
table(8,4)=std(temp) ; 

% compute 25th percentile (first quartile) for ductl 
table(5, 4) = median(temp(find(temp<median(temp)))) ; 

% compute 75th percentile (third quartile) for ductl 
table(7, 4) = median(temp(find(temp>median(temp)))); 

% find mean, median and std. dev. of thickness 
% surface duct 
counter = 1; 

j=i; 

temp=0 ; 
for i=l:X 

if data(counter,15)==1 && data(counter,25)>0 
temp(j)=data(counter,25); 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(9,4)=mean(temp) ; 
table (11,4)=median(temp) ; 
table(13,4)=std(temp) ; 

% compute 25th percentile (first quartile) for surface duct 
table(10, 4) = median(tempi(find(templ<median(temp)))) ; 

% compute 75th percentile (third quartile) for surface duct 
table(12,4) = median(tempi(find(templ>median(temp)))); 
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counter=l ; 

j=i; 

temp=0 ; 
for i=l:X 

if data(counter, 3)==1200 ... . 

&& data{counter,7)==data{counter,8) && data{counter,7)>0 

temp{j)=data{counter, 7) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

a=length(temp); 

table{2,4)=(a/day)*100; % sf duct in day 
counter=l; 

j=i; 

temp=0 ; 
for 1=1:X 

if data(counter,3)==0 .... 

&& data(counter,7)==data(counter,8) && data(counter,7)>0 

temp(j)=data(counter,7); 

j=j+i; 

end 

counter=counter+l ; 

end 

a=length(temp); 

table(3,4)=(a/night)*100; % sf duct in night 

% finds mean, median, first and second quartiles 
% std. deviation of strengt of ducts 
% First duct 
counter = 1; 

j=i; 

for i=l:X 

if data(counter, 32) >0 

tmpl(j)=data(counter, 32 ) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(14,1)=mean(tmpl) ; 
table (16,1)=median(tmpl); 
table(18,l)=std(tmpl) ; 

% compute 25th percentile (first quartile) for ductl 
table(15,1) = median(tmpl(find(tmpl<median(tmpl)))); 

% compute 75th percentile (third quartile) for ductl 
table(17,1) = median(tmpl(find(tmpl>median(tmpl)))); 

% finds mean, median, first and second quartiles 
% std. deviation of strengt of ducts 
% Second duct 
counter = 1; 

j=i; 

for i=l:X 
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if data{counter, 33) >0 

tmp2{j)=data{counter, 33) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table{14,2)=mean{tmp2) ; 
table (16,2)=median{tmp2) ; 
table{18,2)=std{tmp2) ; 

% compute 25th percentile (first quartile) for duct2 
table(15,2) = median{tmp2{find{tmp2<median{tmp2)))); 

% compute 75th percentile (third quartile) for duct2 
table(17,2) = median(tmp2(find(tmp2>median(tmp2)))) ; 

% finds mean, median, first and second quartiles 
% std. deviation of strengt of ducts 
% Third duct 
counter = 1; 

j=i; 

temp3=0 ; 
for i=l:X 

if data(counter, 34) >0 

tmp3(j)=data(counter, 34 ) ; 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(14,3)=mean(tmp3) ; 
table (16,3)=median(tmp3); 
table(18,3)=std(tmp3) ; 

% compute 25th percentile (first quartile) for duct3 
table(15,3) = median(tmp3(find(tmp3<median(tmp3)))); 

% compute 75th percentile (third quartile) for duct3 
table(17,3) = median(tmp3(find(tmp3>median(tmp3)))) ; 


% finds mean, median, first and second quartiles 
% std. deviation of strengt gradient of ducts 
% First duct 
counter = 1; 

j=i; 

tmpl=0 ; 
for i=l:X 

if data(counter, 35) >0 

tmpl(j)=data(counter, 35) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(19,1)=mean(tmpl) ; 
table (21,l)=median(tmpl); 
table(23,l)=std(tmpl) ; 
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% compute 25th percentile (first quartile) for ductl 
table(20,1) = median(tmpl(find(tmpl<median(tmpl)))) ; 

% compute 75th percentile (third quartile) for ductl 
table(22,1) = median(tmpl(find(tmpl>median(tmpl)))); 

% finds mean, median, first and second quartiles 
% std. deviation of strengt gradient of ducts 
% Second duct 
counter = 1; 

j=i; 

tmp2 = 0 ; 
for i=l:X 

if data(counter, 36) >0 

tmp2(j)=data(counter, 3 6); 

j=j+i; 

end 

counter=counter+l ; 

end 

table(19,2)=mean(tmp2) ; 
table (21,2)=median(tmp2) ; 
table(23,2)=std(tmp2) ; 

% compute 25th percentile (first quartile) for duct2 
table(20,2) = median(tmp2(find(tmp2<median(tmp2)))); 

% compute 75th percentile (third quartile) for duct2 
table(22,2) = median(tmp2(find(tmp2>median(tmp2)))) ; 

% finds mean, median, first and second quartiles 
% std. deviation of strengt gradient of ducts 
% Third duct 
counter = 1; 

j=i; 

tmp3=0 ; 
for i=l:X 

if data(counter, 37) >0 

tmp3(j)=data(counter, 31 ) ; 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(19,3)=mean(tmp3); 
table (21,3)=median(tmp3); 
table(23,3)=std(tmp3) ; 

% compute 25th percentile (first quartile) for duct3 
table(20,3) = median(tmp3(find(tmp3<median(tmp3)))); 

% compute 75th percentile (third quartile) for duct3 
table(22,3) = median(tmp3(find(tmp3>median(tmp3)))) ; 
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% finds mean, median, first and second quartiles 
% std. deviation of strengt of ducts 
% surface ducts 
counter = 1; 

j=i; 

tmpl=0 ; 
for i=l:X 

if data{counter,15)==1 

if data(counter, 32) >0 

tmpl{j)=data{counter, 32 ) ; 

j=j+i; 

end 

end 

counter=counter+l; 


end 

table(14,4)=mean(tmpl) ; 
table (16,4)=median(tmpl); 
table{18,4)=std(tmpl) ; 

% compute 25th percentile (first quartile) for ductl 
table(15, 4) = median(tmpl{find{tmpl<median(tmpl)))); 

% compute 75th percentile (third quartile) for ductl 
table(17,4) = median(tmpl(find(tmpl>median(tmpl)))); 

% finds mean, median, first and second quartiles 
% std. deviation of strength GRADIENT of ducts 
% surface ducts 
counter = 1; 

j=i; 

tmpl=0 ; 
for i=l:X 

if data(counter,15)==1 

if data(counter, 35) >0 

tmpl(j)=data(counter, 35) ; 

j=j+i; 

end 

end 

counter=counter+1 ; 

end 

table(19,4)=mean(tmpl) ; 
table (21,4)=median(tmpl) ; 
table(23,4)=std(tmpl) ; 

% compute 25th percentile (first quartile) for ductl 
table(20,4) = median(tmpl(find(tmpl<median(tmpl)))) ; 

% compute 75th percentile (third quartile) for ductl 
table(22, 4) = median(tmpl(find(tmpl>median(tmpl)))); 

% finds mins and maxs of first duct 
counter = 1; 

j=i; 

tempminl=0 ; 
tempmaxl=0 ; 
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for i=l:X 

if data(counter,8)~=0 

tempminl{j)=data{counter, 8) ; 
tempmaxl{j)=data{counter, 6) ; 

j=j+i; 


end 

counter=counter+l ; 

end 

table(24,1)=min(tempminl) ; 
table(25,1)=max(tempmaxl) ; 

% find mins and maxs of second duct 
counter = 1; 

j=i; 

tempmin2 = 0 ; 
tempmax2 = 0 ; 
for 1=1:X 

if data(counter,11)~=0 

tempmin2(j)=data(counter, 11) ; 
tempmax2(j)=data(counter, 9) ; 

j=j+i; 


end 

counter=counter+l ; 

end 

table(24,2)=min(tempmin2) ; 
table(25,2)=max(tempmax2) ; 

% find mins and maxs of third duct 
counter = 1; 

j=i; 

tempmin3=0 ; 
tempmax3=0 ; 
for 1=1:X 

if data(counter,14)~=0 

tempmin3(j)=data(counter, 14 ) ; 
tempmax3(j)=data(counter, 12) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(24,3)=min(tempmin3) ; 
table(25,3)=max(tempmax3) ; 

% find mins and maxs of surface duct 
counter = 1; 

j=i; 

tempmin=0; 
tempmax=0; 
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for i=l:X 

if data(counter,15)==1 && data{counter,8)~=0 

tempmin{j)=data{counter, 8) ; 
tempmax{j)=data{counter, 6) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(24,4)=min(tempmin) ; 
table(25,4)=max(tempmax) ; 

% finds mean, median, first and second quartiles 
% std. deviation of Mdeficit of 
% surface ducts 
counter = 1; 

j=i; 

tempdeficit=0 ; 
for i=l:X 

if data(counter, 39) > 0 

tempdeficit(j)=data(counter, 3 9); 

j=j+i; 

end 

counter=counter+l ; 

end 

table(26,4)=mean(tempdeficit) ; 
table (28,4)=median(tempdeficit) ; 
table(30,4)=std(tempdeficit) ; 

% compute 25th percentile (first quartile) for ductl 
table(27,4) = 

median(tempdeficit(find(tempdeficit<median(tempdeficit)))); 

% compute 75th percentile (third quartile) for ductl 
table(29,4) = 

median(tempdeficit(find(tempdeficit>median(tempdeficit)))); 

%Finds mean, median, lower and higher quartile and std. deviation of 
%surface pressure if there is a duct 
counter = 1; 

j=i; 

temp=0 ; 
for i=l:X 

if data(counter, 6) > 0 

temp(j)=data(counter, 22) ; 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(31,1)=mean(temp); 
table (33,1)=median(temp) ; 
table(35,l)=std(temp) ; 
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% compute 25th percentile (first quartile) for ductl 
table(32,1) = median(temp(find(temp<median(temp)))); 

% compute 75th percentile (third quartile) for ductl 
table(34,1) = median(temp(find(temp>median(temp)))) ; 

%Finds mean, median, lower and higher quartile and std. deviation of 
%surface pressure if there is no duct 
counter = 1; 

j=i; 

temp=0 ; 
for i=l:X 

if data(counter,6)== 0 

temp(j)=data(counter, 22) ; 

j=j+i; 

end 

counter=counter+1 ; 

end 

table(31,2)=mean(temp) ; 
table (33,2)=median(temp) ; 
table(35,2)=std(temp) ; 

% compute 25th percentile (first quartile) for ductl 
table(32,2) = median(temp(find(temp<median(temp)))); 

% compute 75th percentile (third quartile) for ductl 
table(34,2) = median(temp(find(temp>median(temp)))) ; 

%Finds mean, median, lower and higher quartile and std. deviation of 
%surface pressure if there is a surfaceduct 
counter = 1; 

j=i; 

temp=0 ; 
for i=l:X 

if data(counter,15)==1 

temp(j)=data(counter, 22) ; 

j=j+i; 

end 

counter=counter+l ; 

end 

table(31,3)=mean(temp); 
table (33,3)=median(temp) ; 
table(35,3)=std(temp) ; 

% compute 25th percentile (first quartile) for ductl 
table(32,3) = median(temp(find(temp<median(temp)))); 

% compute 75th percentile (third quartile) for ductl 
table(34,3) = median(temp(find(temp>median(temp)))) ; 

%//end of file 


67 


E. 


MATLAB CODE TO FIND DUCTS MONTH-BY-MONTH 


%finds ducts for january 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
monthl=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 
tempcpr= temp==' 01 ' ; 

if tempcpr{2)==1 && tempcpr(3)==1 

monthl{countermonth,:)= data(counter,:); 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1; 

end 

%finds ducts for february 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
month2=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 

tempcpr= temp==' 02 

if tempcpr{2)==1 && tempcpr{3)==1 

month2{countermonth, :)= data(counter, :) ; 
countermonth = countermonth + 1; 

end 

counter = counter + 1 ; 

end 

%finds ducts for march 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
month3=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 

tempcpr= temp==' 03 

if tempcpr{2)==1 && tempcpr{3)==1 

month3{countermonth, :)= data(counter, :) ; 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1 ; 

end 

%finds ducts for april 
x=length(data{: , 1)) ; 
counter = 1; 
countermonth =1; 
month4=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 
tempcpr= temp==' 04 ' ; 

if tempcpr(2)==1 && tempcpr{3)==1 

month4(countermonth, :)= data(counter, : ) ; 
countermonth = countermonth + 1 ; 
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end 

counter = counter + 1 ; 

end 

%finds ducts for may 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
month5=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 
tempcpr= temp==' 05 ' ; 

if tempcpr{2)==1 && tempcpr{3)==1 

months{countermonth, :)= data(counter, :) ; 
countermonth = countermonth + 1; 

end 

counter = counter + 1 ; 

end 

%finds ducts for june 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
month6=[]; 
for i=l:X 

temp = num2str(data(counter,2)) ; 

tempcpr= temp==' 06 

if tempcpr{2)==1 && tempcpr{3)==1 

months{countermonth, :)= data(counter, :) ; 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1 ; 

end 

%finds ducts for july 
x=length(data{:,1)) ; 
counter = 1; 
countermonth =1; 
month7=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 
tempcpr= temp==' 07 ' ; 

if tempcpr(2)==1 && tempcpr(3)==1 

month7{countermonth, :)= data(counter, :) ; 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1 ; 

end 

%finds ducts for august 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
month8=[]; 
for i=l:X 

temp = num2str(data(counter,2)) ; 

tempcpr= temp==' 08 

if tempcpr{2)==1 && tempcpr{3)==1 
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months{countermonth, :)= data(counter, : ) ; 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1; 

end 

%finds ducts for September 
x=length (data{: , 1) ) ; 
counter = 1; 
countermonth =1; 
month9=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 

tempcpr= temp==' 09 

if tempcpr{2)==1 && tempcpr{3)==1 

month9{countermonth, :)= data(counter, : ) ; 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1 ; 

end 

%finds ducts for October 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
monthl0=[]; 
for i=l:X 

temp = num2str(data(counter, 2)) ; 

tempcpr= temp==' 10 

if tempcpr{2)==1 && tempcpr{3)==1 

monthlO{countermonth,:)= data(counter,:); 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1 ; 

end 

%finds ducts for november 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
monthll= [ ]; 
for i=l: X 

temp = num2str(data(counter,2)); 
tempcpr= temp==' 11 ' ; 

if tempcpr{2)==1 && tempcpr{3)==1 

monthll{countermonth,:)= data(counter,:); 
countermonth = countermonth + 1 ; 

end 

counter = counter + 1 ; 

end 

%finds ducts for december 
x=length(data{:,1)); 
counter = 1; 
countermonth =1; 
monthl2= [ ]; 
for 1=1: X 

temp = num2str(data(counter,2)); 
tempcpr= temp==' 12 
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if tempcpr{2)==1 && tempcpr{3)==1 

monthl2{countermonth,:)= data{counter, 
countermonth = countermonth + 1; 

end 

counter = counter + 1 ; 

end 

%//end of file 
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